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Abstract 

The paper provides an introductory discussion about two fundamental models of oscil- 
lator synchronization: the (continuous-time) diffusive model, that dominates the mathe- 
matical literature on synchronization, and the (hybrid) kick model, that accounts for most 
popular examples of synchronization, but for which only few theoretical results exist. The 
paper stresses fundamental differences between the two models, such as the different con- 
traction measures underlying the analysis, as well as important analogies that can be drawn 
' in the limit of weak coupling. 
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Synchronization is a pervasive concept in science and engineering. Currently, it is perhaps the 
most widely studied dynamical concept across systems biology [24,26,83], neuroscience [321135], 
chemistry [39] , physics [3"4"ll62] , astronomy [6j , and engineering [551175] . Because synchronization 
involves interconnection at its core, the relevance of systems theory to model, understand, and 
control synchronization is obvious and was recognized early, e.g. [52J. 

Two fundamental mathematical models of synchronization have emerged across the litera- 
^ . ture: the diffusive model and the kick model (a nickname throughout the paper for pulse-coupled 

synchronization model) . The diffusive model analyzes synchronization as the result of diffusive 
coupling: the interconnection has the input-output interpretation of a static diffusive passive 
map. Owing to the fundamental homogenization nature of diffusion, diffusive interconnections 
tend to reduce differences between the time-course of interconnected variables, thereby favoring 
synchronized behavior if they are strong enough. In contrast, the kick model analyzes syn- 
chronization as the result of mutual rhythmic locking by short and weak pulses, akin to the 
physical phenomenon of resonance. The impulsive nature of the coupling combined with the 
continuous-time flow of the model between the pulses results in a hybrid model, see [51] for a 
rigorous description of the kick model as a hybrid model. 

The diffusive model is largely dominant in the mathematical literature of synchronization. 
Synchronization between trajectories of state-space models is analyzed as an incremental sta- 
bility property [43] (i.e. the trajectories converge to one another rather than being attracted 
toward some equilibrium position) . The leading concepts of Lyapunov analysis [3J , dissipativity 
analysis [1E2E3ESEI] , and— to a growing extent— contraction analysis [MlfMIMIHHllrMIHn] . 
provide natural system theoretic tools to study synchronization. The literature of synchroniza- 
tion (closely related to consensus theory and coordination theory) is growing and the topic has 
attracted many systems and control researchers in the recent years. 

The kick model is largely dominant in natural manifestations of synchronization. Popular 
examples include synchronization of metronomes |84j . clocks [5j[33], heart beats [61] . flashing 
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fireflies [TT], neurons [23], earthquakes |56j . and in fact most if not all spiking oscillators. In 
addition, kick synchronization is a source of inspiration for engineering applications (e.g. syn- 
chronization in wireless sensor networks [31], unsupervised classification problems [B3]). Despite 
the widespread occurrence of the phenomenon, the mathematical literature on kick synchroniza- 
tion is rather sparse compared to the literature on diffusive synchronization. 

Primarily motivated by the recent thesis [44] . the present paper aims at comparing and con- 
trasting the diffusive model and the kick model for the synchronization of periodic oscillators. 
We stress both the differences and the analogies between the two models, with a particular 
emphasis on their global stability properties. The discussion is tutorial in nature and focuses 
on simple examples, such as the coupling of van der Pol oscillators, which provides an insight- 
ful illustration of diffusive synchronization in the weakly nonlinear oscillation regime and of 
kick synchronization in the relaxation oscillation regime. The diffusive model is studied in 
continuous-time models while the kick model, hybrid in nature, is typically studied in discrete 
time. Ultimately, synchronization is always proven by showing that a certain distance between 
trajectories contracts over time, but the contraction measure is distinctively different in diffusive 
and kick models. 

While the diffusive and kick models of oscillator synchronization are fundamentally differ- 
ent, they also exhibit a remarkable analogy in the limit of weak coupling. This is because 
arbitrary oscillator models all reduce to one-dimensional phase models when the interconnec- 
tion is sufficiently weak to maintain system trajectories in the neighborhood of the limit cycle 
oscillations of the uncoupled oscillators. Rooted in the seminal contributions of Winfree [82j and 
Kuramoto [38], phase models of interconnected oscillators have a universal structure entirely 
characterized by their coupling function, which is strongly related to the phase response curve 
of the oscillators (i.e. a function which corresponds to the phase sensitivity of the uncoupled 
oscillators to an external perturbation). As a result, the fundamental difference between dif- 
fusive synchronization and kick synchronization is entirely coded in the shape of the coupling 
function, a phase map defined on the nonlinear unit circle. It is typically harmonic in the weak 
coupling limit of diffusive synchronization and typically monotone (and hence discontinuous) in 
the weak coupling limit of kick synchronization. Again, the synchronization mechanisms and 
the contraction measure are distinctively different even in the weak coupling limit, despite the 
shared model structure. 

The paper structure is illustrated in Figure [TJ Section [2] reviews state-space models of os- 
cillators and their phase reduction. The next sections present global stability results for the 
different synchronization models. Section [3J focuses on (possibly strong) diffusive synchroniza- 
tion while Section 0] focuses (possibly strong) kick synchronization. Section [5] deals with the 
phase models encountered in the limit of weak coupling. Section [S] provides concluding remarks. 

2 Open oscillator models 

This section provides a short introduction to oscillators viewed as open dynamical systems, that 
is, as dynamical systems that interact with their environment [UTj. We first recall basic defini- 
tions about stable periodic orbits in n-dimensional state-space models (see [2TJ[28] for details). 
We then introduce (finite and infinitesimal) phase response curves as fundamental mathematical 
information required for the reduction. We finally show how to reduce n-dimensional state-space 
models into one-dimensional phase models depending on the nature of the input. 
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Figure 1: The paper is organized according to the coupling models. Sections [3] and H] focus 
on (possibly strong) diffusive and kick synchronization, respectively. Section [5] deals with the 
phase models encountered in the limit of weak coupling. 

2.1 State-space models 

We consider open dynamical systems described by nonlinear time-invariant state-space models 

x = F(x) + G(x)u, xel B ,«El, (la) 

y = H(x), y£R, (lb) 

where the vector fields F and G, and the measurement map H support all usual smoothness 
conditions that are necessary for existence and uniqueness of solutions. We write xo, u) for 
the solution of the initial value problem (jlap with x(0) = Xo- 

An oscillator is an open dynamical system whose zero-input steady-state behavior is peri- 
odic rather than constant. Formally, we assume that the zero-input system x = F(x) admits a 
(locally hyperbolic) stable periodic orbit 7 with period T (and angular frequency u) = 2n/T). 
Picking an initial condition Xq on the periodic orbit 7, this latter is described by the (non- 
constant) periodic trajectory <!>(£, Xq,0) = x 7 (t), such that x 7 (t) = x 7 (i + T). The basin of 
attraction of 7 is the maximal open set #(7) from which the periodic orbit attracts. 

Since the periodic orbit 7 is homeomorphic to the unit circle S 1 , it is naturally parametrized 
by a single scalar phase. Any point p 6 7 is associated with a phase 6 € S , such that 

p = x 7 (#M 

(where Xq is by convention associated with the phase 9 = 0). 

For hyperbolic periodic orbits, the notion of phase is extended to any point q in the basin 
of attraction #(7) through the concept of asymptotic phase. The asymptotic phase map : 
#(7) —s- S 1 assigns to each point q in the basin 0(7) its asymptotic phase 9 € S , such that 

lim ||*(i,q,0) - $(t,x 7 (#/w),0)|| 2 = 0. (2) 

t— s-+oo 

This mapping is constructed such that the image of Xq is equal to and such that the progression 
along any orbit in $(7) (in absence of perturbation) produces a constant increase in 9, that is, 
4e(*(t,Xo,0))=a;. 
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Figure 2: The asymptotic phase map : 6(7) — > S 1 assigns to each point q in the basin 6(7) 
a single scalar phase 9 on the unit circle § , such that limt_» +00 ||<fr(i, q, 0) — <&(i, p, 0)|| 2 = 
where p = x 7 (#/u>). The set of all points q characterized by the same phase 6 is the isochron Tq. 
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Figure 3: The van der Pol oscillator exhibits two different oscillation regimes: the quasi- 
harmonic (/i <C 1) and the relaxation (/i 3> 1) oscillation regimes. Quasi-harmonic and re- 
laxation regimes are displayed in (x, x) and (x, z) state-spaces, respectively (with the transfor- 
mation z = x — x 3 /3 — x/n). 

An isochron is a level set of the asymptotic phase map 0, that is, the set of all points in the 
basin of attraction of 7 characterized by a same asymptotic phase. 

The van der Pol oscillator: an illustrative model 

In this paper, we illustrate most concepts on the van der Pol oscillator 

x — — x 2 )x + x = eu, iel, (3) 

where the parameter [i > measures the nonlinearity of the oscillator and the constant e > 
measures the input strength (with |n| < 1 for all times). Historically, this equation modeled 
a simple electrical circuit with nonlinear resistance and was used by van der Pol to study 
oscillations in vacuum tube circuits |77j . It played a seminal role in the development of nonlinear 
oscillation theory. One reason of this success is its ability (with only one parameter) to exhibit 
two very different regimes of oscillations (Figure [3]). For weak nonlinearities (// 1), the 
oscillator displays quasi-harmonic oscillations. For strong nonlinearities (/i ^> 1), it displays 
relaxation oscillations. 

Quasi-harmonic limit 

To study the van der Pol oscillator in the quasi-harmonic limit, it is convenient to rewrite §3§ 
in polar coordinates (x,x) = (r sin(</>), r cos(</>)) as 

r = fj,g(r sin(0), r cos{<p)) cos(</>) + e cos((p)u (4a) 

</>=! — — g(r sm(<f>), r cos(0)) sin(0) sin(^)n (4b) 
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Figure 4: The shape of the infinitesimal phase response curve (for the van der Pol oscillator) 
is very different in both regimes. Typically, it is harmonic in the weakly nonlinear oscillation 
regime and monotone (and hence discontinuous) in the relaxation regime. 



where we denote g(x,x) = (1 — x 2 )x to simplify notations. 

For small values of (i and e (/U,e <C 1), standard averaging theory guarantees that r stays 
in a 0(fi, e)-neighborhood of r* = 2 (see [36] for details). Substituting r by f = 2 + 0(fi,e) 
into (|4bp and keeping first-order terms in the equation yield 

4> ~ 1 — ^g(f sin(f/>), f cos(</>)) sin((/)) — - sin((/))u. (5) 

The approximately equal sign (~) means that (JSJ neglects higher order terms in (/i,e). This 
one-dimensional equation describes the dynamics of the angular coordinate <f> € S 1 . Note that 
this (geometrical) angular coordinate is different from the (temporal) asymptotic phase defined 
in ©. 

Since the angular coordinate dynamics © are one-dimensional, the asymptotic phase map 
appears as a bijective change of variable 6 = 6(</>) given by 

@W : + ~ U I l-^(rsin(6!fcos(e))sin ( e)^ 

This change of variable rescales the state-space and the (temporal) phase dynamics are given 
by 

Q^u + e -ujsmjcj)) 

f — ng{f sin(</>), r cos(0)) sin(</>) 

S v ' 

= :ZQH(0) 

where (j> = Q~ 1 (6). The phase dynamics © are the addition of two terms: the first term 
represents the autonomous angular frequency and the second term represents the influence 
of the input on the dynamics. The function Zqh(-) captures the sensitivity of the oscillator 
phase dynamics to the input. It is known as the (input) infinitesimal phase response curve. 
(This notion will be defined properly in the next section.) For values of fi tending to 0, the 
asymptotic phase map tends to the identity, the angular frequency oj tends to 1, and the 
(input) infinitesimal phase response curve Z QH (e) tends to -isin((9) (Figureg]). 

Relaxation limit (and integrate-and-fire oscillators) 

To study the van der Pol oscillator in the relaxation limit, it is convenient to rewrite ([3]) in 
Lienard's coordinates ( x, z^j as 

— xx' = x — x 3 /3 — z (7a) 
z' = x — eu (7b) 



5 



where we use the transformation z = x — x 3 /3 — x/f-i and where (•)' denotes the derivative with 
respect to s = t/fi. 

For large values of \x (1//-L 2 <C 1), standard singular perturbation theory reduces the dynam- 
ics ([7]) to (see [36] for details) 

x [ 

x' = - 2 - e- (8) 

1 — x z 1 — x z 

on the critical manifold defined by z = x — x 3 /3 and to instantaneous 'jumps' at the folds in the 
critical manifold. Exploiting the central symmetry of the drift vector field (invariance under 
point reflection through the origin), we reduce the dynamics to the one-dimensional dynamics 
on the left branch of the critical manifold: the state x monotonically increases on [—2, —1] 
according to (|8|) and is reset to the lower threshold x = — 2 when reaching the upper threshold 
x = -1. 

Again, since the state dynamics are one-dimensional, there is a bijective change of variable 
9 = 6(x) given by 

i-e 



6(x) :x^u — d£ . (9) 
Jx ? 

This change of variable rescales in such a way that the lower threshold x = —2 is mapped to 
9 = and the upper threshold x = — 1 to 9 = 2ir. The (temporal) phase dynamics are then 
given by 

9' = U } + e(-^)u (10) 

where x = G~ 1 (9). Here again, the phase dynamics ()10p are given by the addition of two terms: 
the autonomous angular frequency and the coupling term. In this case, the phase sensitivity 
function (or infinitesimal phase response curve) Zr(-) is monotone on [0, 2ir) (Figure HJ. 

In the relaxation limit, the van der Pol oscillator is equivalent to an integrate-and-fire model. 
More generally, the integrate-and-fire dynamics are expressed as one-dimensional state dynamics 
between two threshold values (see [HE?]): a scalar state variable x monotonically increases 
between two thresholds x and x, according to the dynamics 

x = F(x), with F(x) > 0, 

for all x G [x, x]. Upon reaching the upper threshold x, the state is instantaneously reset to the 
lower threshold x. Roughly speaking, the oscillator integrates between the two thresholds and 
fires when reaching the upper threshold. 

The most popular integrate-and-fire oscillator is the leaky integrate-and-fire (LIF) oscillator, 
characterized by the monotone vector field F(x) = S+Rx > 0, Vx € [x, x] = [0, 1]. An important 
generalization of the LIF oscillator — in the sense that the dynamics are not monotone anymore — 
is the quadratic integrate-and-fire (QIF) oscillator, defined by the vector field i^x) = S + x 2 , 
with S > [IS]. 

Similarly to ([9]), the asymptotic phase map that corresponds to the integrate-and-fire 
dynamics is the bijective change of variable given by 

f x 1 

G(x) : x i— y uj / , 



_ no 

with the lower threshold x (resp. the upper threshold x) being mapped to 9 = (resp. 9 = 2ir). 

It is worth mentioning that, in the relaxation limit, the van der Pol oscillator model closely 
resemble the popular model of FitzHugh-Nagumo [22|.l51j . a two-dimensional qualitative reduc- 
tion of Hodgkin-Huxley model of neuronal action potentials [30] . Integrate-and-fire models are 
broadly used in neurodynamics [HI37]. 
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2.2 Phase response curves 

For many oscillators, the structure of the asymptotic phase map and therefore the topology of 
isochrons are very complex. This often makes their analytical computation impossible and even 
their numerical computation intractable (at least very expensive for high dimensional oscillator 
models), an issue that prevents from building an exact one-dimensional phase model valid in the 
whole basin of attraction. However, in many situations, a complete knowledge of the isochrons 
is not required to study the oscillator dynamics. Instead, it is sufficient to consider the phase 
response curve, as it has naturally appeared in ([6]) and (fTU|) through the reduction of the van 
der Pol dynamics. 

Starting with the pioneering work of Winfree [82.83J, the phase response curve of an oscillator 
has proven a useful input-output tool to study oscillator dynamics. It indicates how the timing 
of inputs affects the timing (steady-state phase shift) of oscillators. Phase response curves are 
directly related to isochrons but capture only partial information about them. 

Definition 1 The finite Phase Response Curve (PRC) corresponding to a Dirac delta input 
u(-) = e<5(-) is the map Z t : S 1 — > (— ir, ir] defined as 

Z e (9) = Q(^(0/uj) + eG(^ (9 /oj)))- 9(x 7 (6>/w)) . 

after the impulsion before the impulsion 

It associates with each point on periodic orbit (parametrized by its phase 6) the phase shift 
induced by the input. j 

In many situations, the PRC can be determined experimentally. Moreover, it can be obtained 
numerically by computing the perturbed and unperturbed trajectories of the nonlinear state- 
space model and by comparing the asymptotic phase difference between each pair of trajectories. 

A mathematically more abstract — yet very useful — tool is the infinitesimal phase response 
curve, which appears in © and (|10p . It records essentially the same information as the finite 
phase response curve but for infinitesimally small Dirac delta input (e <C 1). 

Definition 2 The (input) infinitesimal Phase Response Curve (iPRC) is the map Z : S 1 — > IR 
defined as the directional derivative 

Z(6) = £>e(xT(0/w))[G(x?(0/a;))] 

where 

„ n , G(x + er?) - 0(x) 
£>6(x) to = lim — — ■ 

e->0 e 

The directional derivative can be computed as the inner product 

L»9(x)[G(x)] = (V x e(x),G(x)) 

where V x O(x) is the gradient of at x and is known as the state infinitesimal phase response 
curve. j 

Remark 1 For small values of e (e<^l), the finite phase response curve is well approximated 
by the infinitesimal phase response curve, that is Z e (-) ~ tZ(-). j 

Remark 2 (Integrate-and-fire oscillators) For integrate- and- fire oscillators, the iPRC has 
the exact analytical expression (see \l(A\3b^ ) 
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Figure 5: For integrate-and-fire oscillators, the finite PRC Z e is directly derived from the 
iPRC Z. 

Moreover, the finite PRC is directly obtained from the iPRC due to the unidimensional nature 
of the state-space. Namely, it follows from Figure [5| that 

Z e (9) = / Z{Q{i))di. 

Jx~<(e/uj) 

For LIF oscillators, the iPRC and the finite PRC are monotone, as for the relaxation van der 
Pol model (Figure^), j 



2.3 Reduced phase models 

We review two popular phase models, which are obtained through phase reduction methods in 
the case of weak input and impulsive input, respectively {TO. 32 , 35 , 39|H2] . 

2.3.1 Weak input 

In the weak perturbation limit, that is, for small input 

u(t) = eu(t), e«l, \u(i)\ < 1 for alii, 

any solution <&(i, xo, u) of the oscillator model which starts in the neighborhood of the hyperbolic 
stable periodic orbit 7 stays in its neighborhood. The n-dimensional state-space model can then 
be approximated by a one-dimensional continuous-time phase model 

9 = lo + eZ{9)u(t) (12) 

where the phase variable 9 evolves on the unit circle S . The phase model is fully characterized 
by the angular frequency u) > and by the iPRC Z : S 1 — > M. 

2.3.2 Impulsive input (kick) 

In the impulsive perturbation limit, the input corresponds to delta-like kicks of amplitude e 
(not necessarily small), that is, 

00 

u(t) = ej^<y(t-t fc ). 

k=0 

Any solution <&(£, xo, u) of the oscillator model which starts from the periodic orbit 7 leaves the 
periodic orbit under the effect of a kick and converges back to the periodic orbit. If the periodic 
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orbit is sufficiently strongly attractive, the trajectory will be back in the neighborhood of the 
periodic orbit before the next kick takes place. The n-dimensional state-space model can then 
be approximated by a one-dimensional hybrid phase model, with 

1. the (constant-time) flow rule 

9 = lo, for all t / t&, (13a) 



2. and the (discrete-time) jump rule (i.e. the kick) 

9 + = 9 + Z e {9), foralU = i fc , (13b) 

where the phase variable 9 evolves on the unit circle S 1 . The phase model is fully characterized 
by the angular frequency co > and by the PRC Z e : S 1 — > (— vr, it]. 

3 Diffusive synchronization 

3.1 Connecting two van der Pol oscillators with a resistor 

Diffusive synchronization is a model of physical interconnection through a diffusive medium. 
As a simple illustration of diffusive synchronization, we consider two van der Pol oscillators 
interconnected with a resistor (Figure E]) 

Xi = -Wi + /j,(xi - xf/3) + Ui (14a) 
Wi = Xi (14b) 
Vi = Xi (14c) 

where Xi and Wi denote the voltage across the capacitor and the current through the inductor, 
respectively. 

The interconnection with a resistor induces a current flow proportional to the voltage dif- 
ference (2/2 — Hi) an d inversely proportional to the resistance R. The smaller the resistance R, 
the higher the coupling strength K = 1/R. Using the vector notation u = (ui,U2) T and 
y = (yi)2/2) T ) the interconnection is expressed as 

u = -Ly 

with the coupling matrix given by 



L = K 



1 -1 
-1 1 



(15) 



This type of coupling is known as diffusive coupling, owing to the nature of the resistor. 

In this context, synchronization is a convergence property for the difference between the 
solutions of different systems. Suppose that we have a network of iV oscillators. The oscillators 
are said to output synchronize if 

lim \\Vi{t)-Vj{t)\\=Q, Vi,j = l,...,N 

t— >+oo 

where || • || denotes the Euclidean norm of the enclosed signal. 

In the following, we show under which conditions the synchronization is guaranteed using: 
incremental stability theory, incremental passivity theory, and contraction theory. 
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Figure 6: The interconnection of two van der Pol oscillators with a resistor causes a current 
flow proportional to the voltage difference (y 2 — y\ ) ■ This interconnection is known as diffusive 
coupling, a name which comes from the diffusive nature of the resistor. 



3.1.1 Incremental stability theory 

Convergence properties for the difference between solutions of a closed system are characterized 
by notions of incremental stability [3] . 
Considering the error variable g — {&xi 

) 1 = xi — x 2 , the error system is written as follows 



2K + n -1 




~A(j)(x) 


1 


e — 






where </>(•) stands for the monotonic function <fi(s) = /xs 3 /3 and A0(-) is defined as A</>(s) = 
<f>(si) — 0(s2). The monotonicity property of the nonlinearity </>(•) implies that 

[0(si) - 0(s 2 )] (si - s 2 ) = A0(s) As > 

for all As = si — s 2 - Considering the Lyapunov function 

V = l(e 2 x + e 2 y )>0, 

we obtain 

V = {n - 2K) e 2 x - A<t>(x) e x < {p. - 2K) e\. 

>o 

Then by the Lyapunov stability theorem and the invariance principle |36t Theorem 4.1 and 4.4], 
we conclude to the synchronization of both oscillators for 2K > fj,. 



3.1.2 Incremental passivity theory 



For open systems, the notion corresponding to incremental stability is incremental dissipativ- 
ity |71j . Denoting the incremental variables by Ax = (Ax x , Ax w ) T = xi — x 2 , Au = u\ — ti 2 , 
and Ay = y\ — y 2 , a system is incrementally passive if it satisfies a dissipative inequality 

AS < w(Au,Ay) 

for an incremental scalar storage function AS 1 (Ax) > with a supply rate w(Au, Ay). 
The incremental system of (|14p is written as follows 



Ax 



Ax + 



Au 



x) 



Ay = [l 0] Ax. 



10 



Considering the incremental storage 

AS(Ax) = ~(A^ + A4)>0, 

we have 

AS = n Ay 2 - Ay A<j>(y) +Au Ay < fj, Ay 2 + Au Ay. 

V v ' 

>o 

Substituting Au = —2KAy in the previous equation yields 

AS < (fi-2K) Ay 2 , 

which implies asymptotic convergence of Ay to zero (that is, output synchronization) when 
2K > «. 



3.1.3 Contraction theory 

Nonlinear contraction theory gives a simple yet general method to study synchronization [80J. 
If the dynamics equations verify 

xi - h(xi ,t) =± 2 - h(x 2 , t) 

where the function h is contracting, then xi and x 2 will converge to each other exponentially. 
Considering the following vector field 

, , ^ \-w + u (x - x 3 /3) - 2K x 
h ( x >*) = ' 

the Jacobian matrix is given by 

Ufjt- 2K) - fix 2 -1 
I 1 

and is negative semidefinite for 2K > u. This implies that h is contracting and that both 
oscillators synchronize when 2K > \i. 



3.2 Large networks 

Most collective phenomena among oscillators in nature arise in large networks of oscillators. 
The notion of diffusive coupling and the tools described in the previous section can be extended 
to a network of N oscillators. 

Each oscillator dynamics is written, for i = 1, . . . , N, as 

Xi = F(xi) + G( Xi )ui, (17a) 
y i = H(x i ). (17b) 

The general diffusive interconnection is then given by 

m = K dvj - y<). * = i, ■ ■ ■ , JV (is) 

where Ku is a positive constant and Ni C J\f is the subset of oscillators transmitting their 
outputs to the ith oscillator. (The set M = {1, . . . , N} denotes all oscillators in the network.) 
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Figure 7: A. Impulsive coupling with integrate-and-fire oscillators: an oscillator reaching the 
upper threshold x = x triggers an instantaneous increment e to the state of the other oscilla- 
tors. B. Impulsive coupling with phase oscillators: an oscillators reaching 8 = 2tt triggers an 
instantaneous increment Z e (9) to the phase of the other oscillators. 



expressed as u = — Ly with the coupling matrix L defined as the Laplacian of the network 
graph 



The diffusive coupling u = — Ly is a passive operator. Neglecting the symmetry neutral 
mode LI = 0, the excess of passivity of the operator is given by the smallest nonzero eigenvalue 
of the symmetric part of L, which is the parameter 2K in (|15p . 

Synchronization is guaranteed in a network of input-output oscillators (fTT|) if the excess of 
passivity in the coupling (2K in the example (|15|) ) compensates for the shortage of incremental 
passivity of the model (|17p (fj, in van der Pol example). A precise statement of this result is 
found in [71] under the assumption of a balanced graph, i.e. when L + L T > 0. The balancing 
assumption has been elegantly removed in the recent paper |13| . 

A limitation of diffusive coupling model is that the incremental stability analysis often 
suggests the necessity of a strong enough coupling. This is in contrast to many synchronization 
problems in which coupling strengths are weak. 

4 Kick synchronization 

4.1 Connecting two van der Pol oscillators with impulses 

A mathematical model of kick synchronization of (identical) oscillators connected through im- 
pulsive coupling was first proposed by Peskin, in the particular case of integrate-and-fire os- 
cillators [61J. When the oscillators fire, they send out a kick that causes an instantaneous 
increment e to the state of all other oscillators of the network (Figure 0A.). 

As a first illustration of Peskin's impulsive coupling, consider two van der Pol oscillators 
in the relaxation limit //> 1. In good approximation, the oscillators are characterized by the 
integrate-and-fire dynamics ©, with x £ [—2, —1]. In addition, suppose that an oscillator which 
fires (i.e. which reaches x = —1) at time sends out a kick u(t) = eS(t — tk) which increases 



Using the vector notations u = (u\ 



. . . ,un) t and y = (yi, . . . ,y^) T , the interconnection is 
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Figure 8: Synchronization of van der Pol oscillators with kick coupling (e > 0). A. The instan- 
taneous state values of an oscillator at the successive firings of the other oscillator approach 
either the lower threshold (x = —2) or the upper threshold (x = —1). (The blue and red sym- 
bols represent the instantaneous state values of the oscillators at the successive firing times.) 
B. After a short transient, the oscillators fire in unison. (The blue and red lines represent the 
firing times of oscillator 1 and 2, respectively.) 
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Figure 9: Phase-locking of van der Pol oscillators with kick coupling (e < 0). A. The instanta- 
neous state values of an oscillator at the successive firings of the other oscillator asymptotically 
converge to a constant value x ~ —1.7. (The blue and red symbols represent the instantaneous 
state values of the oscillators at the successive firing times.) B. The asymptotic firing pattern 
of the oscillators is periodic. (The blue and red lines represent the firing times t& of oscillator 1 
and 2, respectively.) 



the state x of the other oscillator by a value e/(x 2 — 1), according to ©. Through a well- 
chosen change of variable, this increment can be made constant for any state value, so that 
the model is equivalent to Peskin model. Numerical simulations show that the two oscillators 
achieve synchronization: after a short transient period, they fire in unison (Figure EJ. When a 
negative increment e < is considered, the oscillators asymptotically converge to a phase- locked 
configuration: they fire at a constant rate and they are characterized by the same instantaneous 
state values (i.e. x = —2 for the firing oscillator and x ~ —1.7 for the other oscillator) at each 
firing time (Figure [9|) . 

In this paper, we use the more general (but equivalent) definition of Peskin's impulsive 
coupling in terms of phase, a definition which is motivated by the developments of Sections 12.21 
and 12.31 (Figure I7B) . Similarly to (fT3|) . each phase oscillator obeys 

1. the (constant-time) flow rule 

9i = u, if Vj / % : 0j / 2tt, (19a) 
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2. and the (discrete-time) jump rule (i.e. the kick) 



Of = minfJi + Z e (8i), 2vr}, if 3j ^ i : 9, = 2vr, (19b) 

(Note that the coupling is all-to-all, i.e. Mi = M \ {i}-) The threshold imposed in ()19bj) 
corresponds to the absorption phenomenon. If the kick is strong enough, an oscillator may 
aggregate with the oscillator that triggered the kick. The two oscillators have subsequently the 
same phase and create a cluster which behaves as a single oscillator. (We therefore make no 
distinction between a single oscillator and a cluster.) 

Remark 3 The jump rule (|19bj) corresponds to an excitatory coupling, that is Z e (9) > V# £ 
(0, 27r). For the sake of simplicity, we adopt this assumption in the sequel. The extension of the 
results to non- excitatory couplings is straightforward. j 



4.2 Large networks and firing maps 

The mathematical analysis of kick synchronization differs from diffusive synchronization in that 
it can be achieved through the analysis of a discrete-time model. Since the network is uncoupled 
between two kicks, all the information is retained by considering the network state at the discrete 
kick times only. For instance, if the configuration of two coupled oscillators right after a kick is 
(9,2ir), then their configuration right after the next kick is given by (2ir,h(9)), with 

h{9) = 2tt-9 + Z £ (2tt - 9) . 

The discrete-time map 9 + = h(9) expresses the phase differences between the two oscillators at 
the successive kick times. It was originally introduced in |49j as the so-called firing map. 

For the study of large networks, we first assume that the order of the oscillators is not 
modified under the effect of (|19bj) . that is we assume that 9\ + Z e {9\) < 92 + Z e {9i) if 9\ < 92- 
This assumption, which is always satisfied for one-dimensional oscillators such as integrate-and- 
fire oscillators, is summarized as follows: 

Assumption 1 (Order preserving assumption) The finite phase response curve satisfies 
the condition Z' e {8) > -1 V0 G (0,2vr). j 

Provided that Assumption [1] holds, the snapshot configurations of a network of N oscillators 
are given by the successive iterations of a [N — l)-dimensional firing map, which appears as a 
straightforward generalization of the scalar firing map: 



H[(0i 3 . . . , 9n-i)] 



h{8 N „ l -8 l ) 

(20) 

„ h(8N-l — 0JV-2) 



Note that the oscillators are not assigned constant indices but are labeled at each kick according 
to the phase ordering < 8\ < 82 < • • • < 0jV-i < &N = 2vr. 

An extensive study of kick synchronization in large networks is therefore restricted to the 
(global) stability analysis of the firing map (|20p . The strongest stability result is obtained for 
oscillators characterized by a monotone PRC Z e 
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Figure 10: Two features of kick synchronization. A. The oscillators asymptotically converge 
to a unique phase-locked configuration. B. The oscillators achieve full synchronization in finite 
time. 

Theorem 1 Consider a finite PRC that satisfies (i) Assumption^ and (ii) either Z"{9) > 
V0 G (0,2vr) or Z'J(9) < V0 G (0,2vr). Then, the (N - 1)- dimensional firing map $20$, with 
N > 1, has a contraction property with respect to the 1-norm 

N-2 

||(0i,- •• ,0Ar-l)|| = N + ^|0fc-0fc+l| + |0iV-l|- (21) 

k=l 

That is, 

• the firing map is contracting with respect to (|21|) if Z' e {9) < V0 G (0, 2-7r); 

• i/ie firing map is expanding with respect to (|2ip i/ ^(0) > V0 G (0, 2ir). j 

A straightforward corollary of Theorem [1] is that kick synchronization has two important 
features: (i) isolated phase-locked configuration and (ii) finite-time synchronization. 

(i) Isolated phase-locked configuration: If the phase response Z e is monotone decreasing, the 
network globally converges toward the unique fixed point of the (N — l)-dimensional 
firing map, which corresponds to the unique phase-locked configuration of N oscillators 
(Figure fTUR). This behavior was previously obtained in Figure with two van der Pol 
oscillators in the relaxation limit. Except in the weak coupling limit, the phase- locked 
configuration is not a splay state (that is, the phase differences between the successive 
oscillators are not identical). In addition, since a single oscillator may represent a cluster 
of (locally synchronized) oscillators, this configuration also corresponds to a phase-locked 
clustering configuration. 

(ii) Finite-time synchronization: If the phase response Z e is monotone increasing, Theorem [JJ 
implies that the fixed point is (globally) unstable. Then, successive absorptions lead the 
network to full synchronization in finite time, that is, all oscillators share the same phase 
(Figure [TUB) (see also [49] ). This behavior was previously obtained in Figure [8] with two 
van der Pol oscillators in the relaxation limit. 

The study of kick synchronization for general phase dynamics is still (very) limited. To 
our knowledge, there is so far no global result for oscillators that are not characterized by a 
monotone PRC. While some oscillators satisfy the hypotheses of Theorem[JJ(e.g. LIF oscillators, 
van der Pol oscillators in the relaxation limit /i > 1), other oscillators are not characterized 
by a monotone PRC (e.g. QIF oscillators). In this latter case, even a local stability analysis 
may become elusive — although Assumption [1] ensures that the firing map has still a unique 
fixed point — and the networks can display other (more complex) collective behaviors [35]. In 
addition, the general study of oscillators that do not satisfy Assumption [JJ — for which a firing 
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map cannot be denned — remains an open problem. These few examples are all relevant research 
perspectives. 

Whereas few studies have investigated kick synchronization for oscillators with general 
phase dynamics, several extensions of the original impulsive coupling can be found in liter- 
ature (e.g. reduced interconnectivity [TBHTSl [2511581 171 1 176]. delays [201E5], non-instantaneous 
interactions 0E1E21IH5I , non-identical oscillators [7l H2l[I5l lBBj). 

4.3 Infinite populations 

Kick synchronization can also be studied in the continuous limit of TV" — > oo oscillators. In 
this case, a continuum of oscillators is described by a phase density function p(9, t) normalized 
on S 1 , or is equivalently described by a flux 

J(9,t) = p(9,t)v(9,t), 

where v(6, t) is the velocity of the oscillators. The evolution of the oscillators obeys the conti- 
nuity equation 

!Lp{0,t) = ~J{p,t), (22) 

with the boundary condition J(0, t) = J(2tt, t) Vt. 

When the population is infinite, the impulsive coupling is continuous and proportional to 
the flux J(2tt, t). We derive the result as follows. In the case of finite populations, the flow 
rule (|19ap and the jump rule (|19bp imply a velocity 

oo 

v{O i ,t) = u> + Z € iO i ) Y, Z^-'fc )' ( 23 ) 

jeAf\{i} k=0 

(i) 

where tY' denote the times at which oscillator j reaches the phase 6 = 2ir. Since the flux is 

given by J(2ir, t) = 1/N J2jeAf Sfclo ~~ )> ^ f° uows that, for a large number of oscillators 
where N \ {i} ~ M, one has 

v(9i,t)*iu + Z e (0i)NJ(2ir,t). (24) 

In the limit of an infinite number of oscillators, the impulsive coupling is an infinite sum of 
infinitesimal kicks e = K/N <C 1, where K is a positive constant. Then, Remark [1] and (|24p 
imply that 

v(9,t)=uj + KZ(9)J(2iT,t) (25) 

and the coupling is proportional to the (continuous) flux J(2-K,t). 

As a parallel to the results obtained for finite populations, the continuity equation (|22p - (|25p 
has strong stability properties if the iPRC is monotone. Its global stability is shown using the 
continuous analog of the 1-norm (|2ip . which has the interpretation of a total variation distance. 
The result is summarized as follows |47| : 

Theorem 2 Consider an iPRC that satisfies Z"{9) > V0 6 (0, 2?r) or Z"{9) < W € (0, 2vr). 
Then, an admissible solution of (|22p ~ (|25p 

• exponentially converges to the unique stationary solution (when it exists) if Z'(9) < 
V0 G (0,2vr); 

• reaches synchronization (infinite flux) in finite time if Z'{9) > V# G (0, 2tt). j 
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Figure 11: Two collective behaviors for kick synchronization in infinite populations. A. The 
continuum converges to the stationary solution characterized by J (8, t) = J* . B. The continuum 
reaches synchronization in finite time. 

According to Theorem [21 kick synchronization in infinite populations involves two collective 
behaviors, which are the exact analogs of the behaviors observed for finite populations. If the 
iPRC is monotone decreasing, the network converges to the unique stationary solution, which 
corresponds to a constant flux J (9) = J* (Figure fTTTA). This is the analog of the phase- locked 
clustering configuration. If the iPRC is monotone increasing, the flux tends to a Dirac function 
and the network achieves synchronization in finite time (Figure [lib). 

LIF oscillators and van der Pol oscillators in the relaxation limit [i 3> 1 are characterized 
by a monotone iPRC that satisfies the hypotheses of Theorem [2[ If the iPRC is not monotone, 
there is so far no global stability result for the continuity equation (|22p -(|25 p . However, local 
stability of the stationary solution can be studied numerically (see e.g. [75]). 

5 Phase models in the weak coupling limit 

The results presented in the previous sections emphasize important differences between diffusive 
synchronization and impulsive synchronization. Both the analysis techniques and the contrac- 
tion measures differ in a fundamental way. However, it is remarkable that both frameworks are 
unified in the weak coupling limit. Indeed, using averaging techniques, Kuramoto showed that 
the phase dynamics of weakly coupled oscillators can always be reduced to the unique canonical 
form 

e i = u+^2T ij (0 i -e j ), (26) 

where the Tjj(-) are coupling functions closely related to the iPRC [39]. As a consequence, the 
phase dynamics (|26p are a general paradigm that holds both for diffusively coupled oscillators 
and impulsively coupled oscillators. 

5.1 Weak diffusive coupling 

In the case of a weak diffusive coupling, we show that the phase dynamics (|26p correspond to 
the averaged dynamics of the oscillators. (We follow similar lines as in |39| (see also |32|)). 
It follows from (j!2p and (j!8|) that the oscillators are characterized by the phase dynamics 

6\ = u + eZ(6i) KjiiHjiOj) - H^)) , 

with Hi(6i) = Hi{-Ki '{Oi/ui)). The phases can be decomposed as 6i = cot + ipi, where ipi are 
slow phase deviations from the uniform natural oscillation u) t. Then, the phase dynamics are 
rewritten as 

ipi = eZ{ujt + ipi) Y KjiiHjiut + ipj) - Hi(ojt + ipi)) . 
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Note that Z is considered here as the 27r-periodic extension of the iPRC on the real line, 
i.e. Z{x) = Z(x mod 2ir). Next, averaging the above dynamics over a period T = 2it/oj and 
under fixed ipi and i/jj, we obtain 



eK f 

^i=Y\ -tft / Z{ut + ipi) x (Hj(ut + il>j)- Hi{oot + iPi)) dt 

eK f 2n 

y —£ / Ztyi - i>j + a) x (Hj(s) - Hityi - + a)) ds 

7?tr. u JO 



where we have used the change of variable wt + ipj = s. With the coupling functions 

rii(-) = — / Z{- + s) {H^a) - HA + a)) ds , (27) 
w Jo 

the last equation yields 

J^Itf^-^), (28) 

which is equivalent to ([26]) . 



Remark 4 (van der Pol oscillator) Mimicking the computation steps to transform to 
([6]), u>e can write (j!4j) in po/ar coordinates (x,w) = (r sin(</>), — r cos(0)) and appZy averaging 
theory to obtain the phase dynamics given by 

■ u!Cos(4>) 
= lo + e . — - u 

r + fig{r sm(<p)) cos{(p) 

* „ ' 

=:Z(0) 

with (ft given by the appropriate bijective change of variable Q~ 1 (9), f = 2 + 0(fj,,e), and 
g{x) = {x-x 3 /3). 

For values of fi tending to 0, the oscillator output is given by H(6) = 2sin(#) and the 
iPRC by Z{0) = ±cos(0). Applying fl2?J and using basic trigonometry, the coupling function 
is given by rVj (0) = — eKjiir sin(9) . The diffusive interconnection of quasi-harmonic van der 
Pol oscillators leads thus naturally to the popular Kuramoto model characterized by a sinusoidal 
coupling function. j 

5.2 Weak impulsive coupling 

Next, we derive the averaged phase dynamics (|26|) in the case of a weak impulsive coupling (see 
also [H]). From CO) and (J23D, one has 

oo 

9i = u + tZipi) K ji J2 S( t - *fc } ) ' 
jeMi fc=o 

where, in full generality, we have introduced the constants K„ and a general interconnection 
topology J\T%. As in Section l57l] using the phase deviation and averaging the dynamics under 
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a constant ipi yield 



^ = E ^ r z ^ + vo E *(* - *S°) * 

T 



= E ! f i ^-^)> 

where the last equality holds since cut^' + ipj = 2tt and given the periodicity of the iPRC. With 
the coupling functions 

f K 

T ij (-) = ^Z(-), (29) 

we recover (|28p , which corresponds to the general form (|26p • In the case of a (weak) impulsive 
coupling, it is remarkable that the coupling functions are proportional to the iPRC itself. 



5.3 Global stability results and collective behaviors 

In the weak coupling limit, the collective behaviors of diffusively coupled and impulsively coupled 
oscillators depend on the (global) stability properties of the averaged dynamics (|2f)j) . When the 
oscillators are initialized within a semicircle, a change of coordinate maps the dynamics into M. N 
and the analysis of (|26[) is equivalent to a consensus problem on a convex set (see e.g. |50j). 
In this case, a global stability analysis can be performed using consensus theory [50J. In this 
section, we rather discuss the global stability properties of (|26p on the whole (nonconvex) 
torus T , but for the particular all-to-all topology with identical connections (i.e. Mi = J\f \ {i} 
and Kji = K Vi,j). 



5.3.1 Kuramoto model 

The most popular model of the form (|26|) is the Kuramoto model |40j . Characterized by the 
coupling function T{6) = —K/Ns'm(9), which can be considered as the first Fourier harmonic 
of a more complex coupling function, the Kuramoto model appears as the generic model for the 
averaged dynamics of (diffusively) coupled oscillators. 

The Kuramoto model is characterized by the following collective behaviors (for almost every 
initial condition): 

• If K < (inhibitory coupling), the oscillators converge toward an incoherent state char- 
acterized by YlkeAf el6k = (with i = V~ 1) (balanced state); 

• If K > (excitatory coupling), the oscillators achieve full synchronization. 

These global properties result from the fact that (|26[) is a gradient system for Kuramoto 
model [32j . In the case of shifted coupling functions T(9) = sin(# — 5), a general Lyapunov 
function still enforces global convergence toward synchronization or incoherent state |81j . 



5.3.2 Monotone coupling function 

It follows from (|29p that a monotone coupling function is obtained in the case of an impul- 
sive coupling, for oscillators characterized by a monotone iPRC (e.g. LIF oscillators, van der 
Pol oscillators in the relaxation limit /i > 1). (The coupling function is monotone on (0, 2ir) 
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and characterized by a discontinuity r(0 - ) 7^ T(0 + ), since it satisfies the 27r-periodicity condi- 
tion.) Then, as was shown in Section [U the monotonicity property of the oscillators induces a 
contraction property for the averaged dynamics (|26|) . 

The contraction is shown in a rotating frame associated with an oscillator, that is, for the 
equivalent dynamics 

N-l AT-l 

0=1 3=1 

that are obtained by using the change of variable (fi = 9i + \ — 9\ in ([26]) . Note that ([30]) is 
defined in the closure of the cone C = {(ipi,--- , ipN-l) G (0, 2n) N ~ 1 \(pi < tpi+i). The result, 
which is the analog of Theorem [lj is summarized in the following theorem [38] . 

Theorem 3 Consider a monotone coupling function that satisfies T"(9) > V6* G (0,2-7r) or 
T"{9) < V# € (0,2vr). Then, 

• the dynamics (|3Q|) are contracting inC with respect to ()21[) ifT'{9) < V# 6 (0, 2tt); 

• the dynamics (|30p are expanding in C with respect to (|2ip ifT'(9) > V# G (0,27r). j 

A corollary of the contraction property of Theorem [3] is that networks of oscillators coupled 
through a monotone coupling function display two behaviors (for almost every initial condition) : 

• If r is monotone decreasing, the oscillators converge to the unique incoherent configuration 
ip* k = kjf (splay state); 

• If r is monotone increasing, the oscillators achieve full synchronization in finite time. 

This global behavior is similar to the global behavior of Kuramoto oscillators, but also 
characteristic of kick synchronization. The synchronization takes place in finite time, and the 
incoherent state is an isolated fixed point (splay state), whereas it is a (N — 3)-dimensional 
manifold in Kuramoto model. A consequence of this difference is that the asymptotic behavior 
of (j26|) is robust to small heterogeneity in the natural frequencies uj in monotone firing oscilla- 
tors [36] whereas the asymptotic dynamics of Kuramoto model can be highly complex even for 
small heterogeneities |72| . 

5.3.3 Other coupling function 

As for pulse-coupled models, there are only few global stability results for phase models char- 
acterized by generic — non-sinusoidal, non-monotone — coupling functions. When the coupling 
function is odd, the phase oscillators correspond to a gradient system and are characterized by 
global properties similar to Kuramoto model [32]. For more general coupling functions, however, 
a global stability analysis is usually elusive (e.g. weakly pulse-coupled QIF oscillators [33])- In 
this case, a local stability analysis can be performed and, in particular, there exists a criterion 
for the local stability of the incoherent state 0(31] ■ In this context, local stability of clustering 
configurations has been investigated in |27p55j . for instance, and the design of the coupling func- 
tion to achieve particular cluster states is considered in [57] . Note also that several studies have 
considered generic coupling functions in the case of non- identical phase oscillators (e.g. |14|). 
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6 Conclusion 



Diffusive synchronization and kick synchronization are two distinctively different models of syn- 
chronization, underlying different physical synchronization mechanisms and leading to different 
analysis tools. While most popular manifestations of synchronization seem more akin to the 
kick model than to the diffusive model, the literature on kick synchronization is sparse, probably 
owing to the hybrid nature of the kick model and to the mathematical difficulty of analyzing 
nonlinear resonance. For instance, a number of kick synchronization problems remain unsolved 
(e.g. non-identical oscillators) and several research themes are still unexplored (e.g. general — 
application oriented — dynamics) . 

By essence, both the diffusive model and the kick model are crude idealizations of the 
complex synchronization phenomena observed in nature. But it is fair to recognize that the 
importance of the 'kick' deserves more consideration in a mathematical literature dominated 
by the diffusive model. This is for instance illustrated in the recent paper celebrating the 300th 
anniversary of the first scientific investigation of synchronization by Huygens [3]. In that sense, 
the present paper is an invitation to the growing hybrid systems community to contribute — as 
in the recent work |54j — to a deeper understanding of a fundamental property of interconnected 
nonlinear dynamical systems. 
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